Astronomy & Astrophysics manuscript no. 


February 2, 2008 


(DOI: will be inserted by hand later) 





Temporal evolution of magnetic molecular shocks 

I. Moving grid simulations 

p. Lesaffre^'^'^, J.-P. Chieze^, S. Cabrit^, and G. Pineau des Forets^'^ 

Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK 
University of Oxford, Department of Astrophysics, Oxford 0X1 3RH, UK 
CEA/DAPNIA/SAp Orme des Merisiers, F-91191 Gif-sur-Yvette Cedex 

LERMA, UMR 8112 du CNRS, Observatoire de Paris, 61 Av. de I'Observatoire, F-75014 Paris 
IAS, UMR-8617 du CNRS, Universite Paris-Sud, bat. 121, F-91405 Orsay 
LUTH, UMR-8102 du CNRS, Observatoire de Paris, F-92190 Meudon Cedex 



Received September 15, 1996; Accepted March 16, 1997 

Abstract. We present time-dependent ID simulations of multifluid magnetic shocks with chemistry resolved down 
to the mean free path. They are obtained with an adaptive moving grid implemented with an implicit scheme. We 
examine a broad range of parameters relevant to conditions in dense molecular clouds, with preshock densities 
10^ cm"^ < n < lO'' cm ^, velocitie s 10 kms ^ < it < 40 kms ^, and three different scalings for the transverse 
magnetic field : B = 0, 0.1, 1 /iGx -^n/cm^^. 

We first use this study to validate the results of Chieze, Pineau des Forets & Flower (1998), in particular the 
long delays necessary to obtain steady C-type shocks, and we provide evolutionary time-scales for a much greater 
range of parameters. 

We also present the first time-dependent models of dissociative shocks with a magnetic precursor, including the 
first models of stationary CJ shocks in molecular conditions. We find that the maximum speed for steady C-type 
shocks is reached before the occurrence of a sonic point in the neutral fluid, unlike previously thought. As a result, 
the maximum speed for C-shocks is lower than previously believed. 

Finally, we find a large amplitude bouncing instability in J-type fronts near the H2 dissociation limit {u ~ 
25 — 30 kms~^), driven by II2 dissociation/reformation. At higher speeds, we find an oscillatory behaviour of 
short period and small amplitude linked to coUisional ionisation of H. Both instabilities are suppressed after some 
time when a magnetic field is present. 

In a companion paper, we use the present simulations to validate a new semi-analytical construction method for 
young low-velocity magnetic shocks based on truncated steady-state models. 

Key words, magnetohydrodynamics - interstellar matter - shock wave propagation - time dependence - hydrogen 
- molecular oscillations 



1. Introduction 

Today, essentially all observational diagnostics of shocks 
are based on steady state models. Indeed, the steady 
state form of the monodimensional equations of hydro- 
dynamics for these shocks is an ordinary differential 
equation. Therefore, computation of such models can 
encompass the luxury of details involved in the ther- 
mal, magne tic and chemical richness of the interstel- 
lar medium. iHollenbach fc McKeel lll979t 1X98(1 Il989t) de- 
scribe the destruction and reformation of molecules in 
very strong jump (hereafter J-type) sh ocks, after hav- 
ing s t udied th e ir rad iative precursor ijShull fc McKed 
I1979I) . iMiillanI lll97lh discovered that high magnetic 
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fields can transfer kinetic energy to the rmal energy i n 
a continuous mann er (C -type shocks). iDraind lll980l) . 
iDraine et ai\ l|l983l) . and iRoberee fc Draind l)l990|) de- 
scribe the multifluid nature of these shocks. In a se- 
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carefully examine the 
general way, all colli- 


relevant chemistry, and in 
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sional interactions between the charged and neutral flu- 
ids. They point out the strong influence of chemistry on 
the magnetic precursor, carefully address the energy losses 
linked with this chemistry, and stress the need to follow 
the level populations of II2 along the flow. H2 is indeed one 
of the most efficient coolants, and therefore one of the most 
powerful diagnostics. A good review of the steady state J- 
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type a nd C-type models can be found in iDraine fc McKed 
l)l993|) . All these studies assume a purely tra nsverse mag- 
netic field frozen in the ions. Other authors ( Piliop et al. 
1990HPilipp fc Hartauistlll994HCaseni aZ.t.l997;.Wardl& 



1998ft have investigated the dynamical role of the grains 
as well as the effect of oblique magnetic fields. For high 
pre-shock densities, they showed that the drift velocity 
had been underestimated, producing a significant rise in 
the maximum temperature. They also discovered that ro- 
tation of the magnetic field could bring one of the charged 
fluids at rest in the shock frame. 

Time-dependent studies of dense shocks in one dimen- 
sion have focused separately on three points. The first one 
is the oscillatory instabilities due t o the temperature de- 
pend e nce of the coo l ing of t he gas JChevalier fc Irnamural 
1982'; 'Gaetz et al! '1988; 'Strickland & Blondin' 'l995| 
-Wai der & Fo lini ,1996; Li m et aL .2002; Smith .fc Rosen 



The second one is the early evolution of multi- 
fluid shocks, which is found to involve a combi nation 
of C and J shock fronts ( Pineau des Forets et all 

i i( 



1997t ISmith fc Mac Lo^ Il997t IChieze et all Il99. 



Pineau des Forets et al .1 lll997l) and IChieze et al\ l)l99S 



include non-equilibrium cooling and chemistry in the 
fluids, with a degree of details comparable to steady state 
models. They prove that in low ionised magnetised media, 
the steady state could be reached only after very late ages 
(up to a few times 10^ yr), far greater than the variation 
time scales of shock driving sources. The third point is 
the dynamical effects of grains, which is beginning to be 
included in tirn e-dependent studies ijCiolek fc Robergel 
l2002HFallell20f)l . 

Chieze, Pineau des Forets & Flower (1998) used a time 
stretching method to resolve the sharp J-fronts in time- 
dependent magnetised shocks. However, their method did 
not allow to resolve non-viscous fronts (e.g. H2 reforma- 
tion regions), and could introduce synchronisation prob- 
lems. To overcome both limitations, we have developed a 
new moving grid algorithm which allows one to resolve 
discontinuities down to the mean free path. Although we 
do not include the treatment of gra in dynamics to remain 
consistent with lChieze et al\ l)l998|) . we note that our nu- 
merical scheme already provides th e requ i red fr amework 
to solve the problems mentioned bv iFallel l)200,'^ . as it is 
implicit with very high resolution. As such, it is the first 
algorithm able to model time-dependent multifluid disso- 
ciative shocks. 

In this paper, we use our code to validate the evo- 
lutionary behaviour and time scales obtaine d with the 
" anamorphosis" method of lChieze et all l)l998(l , and to ex- 
tend their range of shock parameters to cover the denser 
conditions encountered in protostellar jets and molecular 
clouds, including dissociative and partly ionising shocks. 
We (1) compute the evolution time-scales of these shocks, 
and interpret physically their behaviour, (2) find that sta- 
tionary CJ shocks are obtained before the occurence of 
a sonic point in the shock frame and are thus more fre- 
quent than previously thought, (3) describe two types of 



oscillations driven by chemistry that were not identified 
previously. 

In Sect. 2, we briefly present our num erical model, 
with the improvements made compared to IChieze et all 
l(l998l) . Section 3 describes the results and compares them 
to previous work. Section 4 discusses the limitations of our 
method and Sect. 5 summarises our conclusions. 



2. Numerical and physical inputs 

2.1. Numerical scheme 

2.1.1. Time-step integration 

We use a fully non-linear implicit scheme for time integra- 
tion. The implicitation parameter is set equal to 0.55 ; that 
way, we combine a quasi order 2 accuracy in time with an 
unconditional stability of the scheme. We use a Van-Leer 
advection scheme. The relative variations of the variables 
is kept under 5% at each time-step. We reproduce classi- 
cal tests such as the Sod's shock tube, Rankine-Hugoniot 
relations, and Sedov explosion w ith a 0.3% accu racy. 

More details can be found in iLesaffrel ll20 2 [). 

We solve for the same equations as Ichieze et all 
(1998), although the hydrodynamics and chemistry are 
now solved simultaneously. This amounts to a number of 
33 chemical equations, 2 momentum equations, 4 energy 
equations and one equation for the moving mesh all cou- 
pled together. Solving for the chemistry without splitting 
it from the hydrodynamics helps numerical convergence, 
but has a computational cost, since the Jacobian matrix is 
far heavier to be inverted. In a few cases, the code stalled 
because a proper convergence of the Newton-Raphson it- 
erations required a time step much too small. 

2.1.2. Adaptive grid 

Chemistry encompasses an extensive range of time scales, 
down to times as short as a few mean collision times of the 
particles. When this range of time scales is coupled with 
hydrodynamics, it generates a full range of spatial scales, 
down to a few mean free paths. 

To achieve this extremely high resolution where 
needed, we us e the moving grid algorithm designed by 
IPorfi fc Drurvl l)l987l) . We define the delay time parame- 
ter of the grid as the sound crossing time of the smallest 
zone. Our resolution function combines three arguments ; 
the gradient of the neutral temperature, the gradient of 
the ionic temperature, and the maximum of all chemical 
gradients. 

Thanks to this method, we were able to resolve all 
shocks and chemical fronts with a total number of zones 
as low as 100, wi thout having to res ort to the anamor- 
phosis method of IChieze et all l)l998l) . The present code 
can safely handle two or more sharp features like a strong 
shock followed by a molecular reformation region, or a 
shock, a contact discontinuity and a reverse shock. The 
natural viscosity is used in the shock fronts, resulting in 
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a viscous spread of typically 10 pc/(riH-cm'^) resolved 
with about 10 zones. 



2.1.3. Diffusion 

We consider diffusion for the chemical species, based on a 
Fick law with a diffusion length equal to the local mean 
free path Ad = l/crriH where a = 10~^^ cm^ is the molec- 
ular cross section, and nu is the numerical density of 
hydrogen nuclei. This is motivated by the fact that in 
regions where molecules are reforming, the fluid veloci- 
ties are usually very low. The reformation times of the 
molecules would then lead to an unrealistically thin front. 
This diffusion term ensures that the front has at least a 
width of the order of the mean free path. 



2.2. Physical inputs 
2.2.1. Chemistry 

The network comprises around 120 reactions involving 33 
differen t species (in cluding electrons). It is the same net- 
work as lChieze et all l(l99ii) with a few additional (or up- 
dated) reactions : 

— collisional dissociations of H2 by e~ fFlower et aZ.' 
j"l996) , H, H2, He, and streaming ions (fwilgenbus et ai 

— collisional ionisation of H by e~ I Le Bourlot et a/1 
[2002). 

~ updated H2 formation on grain s : for the sticking coef- 
ficient , we use expressio n (4) of| Hollenbach fc Sa lpeterl 
lll97lh cahbrated by iBuch fc Zhang (,1991.1 and 
iMasuda et~ld\ ((1993). 

To compute the energy transfers and coolings due to 
chemical reactions, we properly track the fluid to which 
each species belong. This is critical especially for endother- 
mic reactions like collisional ionisation of H. 



the same excitation coefhcients for H-CO collisions. 
The optical depth parameter is computed in the static 
approximation, with a distance scale corresponding to 
= 10 mag. 

2.3. Parameter space and initial conditions 

The three parameters n (number density of hydrogen 
nuclei), u (initial velocity), and B (initial ambient 
magnetic field) characterise a simulation. We replace B 
by the parameter b, setting B = &(n.cm^)2 fj,G. Typical 
conditions in the interstellar medium correspond to b = I. 
But as we consider only the transverse component of the 
field, b can take values between and 1. We combine 
n = 10^, 10'*, 10^ cm-^, u = 10, 20, 30, 40 kms"*, 
and 5 = 0, 0.1, 1 and produce a grid of 36 different 
simulations. To test the dynamical instability that we 
observed in dissociative cases, we ran two additional cases 



with 6 = 0, 



25 kms 



lO" and IQ^ cm-3. The 



abundances of He, C, O, and Fe nuclei relative to rin are 
the following : He=0.1 ; C=1.7 10"'' ; 0=4.25 lO"'' ; 
Fe=1.8 lO-'^. 

We use a piston-like protocol. The initial conditions 
of the simulation box are homogeneous, in thermal and 
chemical equilibrium with a density corresponding to a 
given nn = n. The transverse magnetic field is uniform, 
equal to B. All interfaces of the cells in the simulation have 
a velocity u, except the last one, which is a fixed reflexive 
boundary. The moving grid algorithm allows us to begin 
with a very small simulation box of a few mean free paths 
(typically 30). We make the left border of the box contin- 
uously flee the shock front. That way, the adiabatic shock 
is always fully resolved. The simulation is stopped when 
the box is too large and the resolution of sharp fronts can- 
not be supported any more by the algorithm. This occurs 
when the dynamical range of length scales is greater than 
typically a billion. 



2.2.2. Atomic and molecular cooling 

We use refined versi ons of al l the atomic and molecular 
coolings used bv iChieze et al\ (|1998^) . with the addition of 
H cooling. 

At each time-step, we compute the level population 
of C (4 levels), C+ (4 levels) and O (5 levels) excited 
by e~, H, He, a nd H^. Co l lisiona l excitation coefficients 
are taken from Mendoza fl983"); 'Haves fc Nussbaumei' 
1984 ): i Monteiro fc Flower. ( J,987) ; HoUenbach fc McKea 
1989t): iLavallevI ^l200n^ . Lyman a cooling is taken from 
iFlower et al\ |(199(t| ). and cooling due to radiative recom- 
bination of H is from lSnitze r (1978) (case B). We properly 
track the energy lost by each fluid, depending on the col- 
lider. 

We consi der cooling by H2 llLe Bourlot et a/1 
I1999I) OH llHollenbach fc McKej Il979l) . CO and 

H2O ("Neufe ld fc Kaufmanlll99 3^. Since the latter authors 
provide data for excitation of CO by H2 only, we assumed 



3. Results 

3.1. Final steady-states and trajectories in the piston 
frame 

Table ^ summarises the final steady-state of each model 
as a function of the initial parameters b, n, and u. When 
& = 0, all flows evolve to stable J- type shocks, except at in- 
termediate velocities for which we find strong or moderate 
undamped oscillations. When b is finite, the flow generally 
evolves to a steady C-type shock. At sufficiently low b or 
high u, however, a strong J shock in the neutrals persists 
behind the continuous magnetic precursor, and a steady 
CJ-type shock is obtained. Figures 1-2-3 show the thermal 
and chemical structure of typical shocks in their final J, 
CJ and C-type steady-state. Figures 4-5 show intermedi- 
ate states of J-shocks with higher entrance velocities. 

At the bottom of these figures, we also plot the shock 
trajectories, i.e. the positions as a function of time of the 
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C precursor and the J front (defined as the zones with 
maximum ratio of viscous to thermal pressure in the ionic 
and neutral fluids, respectively). We also plot their corre- 
sponding propagation velocity v away from the piston, ob- 
tained by differentiating with respect to time, and smooth- 
ing with a median filter (20 time steps of width). This 
avoids spurious high velocities due to change in the zone 
number. 

Even at intermediate times, v for stable J-type fronts 
is found to obey the steady-state relation v = u/{C — 1) 
where C is the neutral fluid compression factor at the 
piston. In the initial adiabatic phase, v — u/3 {C — i for 
7 = 5/3). It then decreases during the cooling phase, down 
to a steady value when thermal equilibrium is reached in 
the post-shock gas. These three phases can be very clearly 
seen in Fig. Id. For the C-type precursor, the initial veloc- 
ity corresponds to the ion magnetosonic speed, v then de- 
creases as u/ (C —1) where C is now t he magnetic compres- 
sion factor (see iLesaffre et ai\\2004. hereafter Paper II). 
At early times, the J-type feature of a magnetised shock 
behaves like the same shock without magnetic field, as 
neutrals and ions are still decoupled. Thus, the J-front 
propagates much slower than the C-front, which allows the 
magnetic precursor to grow. In the steady C-type case, the 
J-front velocity eventually catches up with the C-front due 
to the softening of its entrance conditions. In the steady 
CJ-type case, the C-front velocity is forced to slow down 
to the J-front value by early effective recoupling between 
the neutrals and the magnetic field through their coUi- 
sional interactions with ions, and the precursor growth is 
halted. 

3.2. Time scales and length scales 

Table m gives, next to the final state of the shock (J, CJ, 
or C), the time scales in years taken to reach steady state. 
It includes in parentheses the width (in pc) of the steady- 
state shock structure (we give the length of both C and J 
components for CJ-type shocks). 

We compared a time evolution sequence for param- 
et ers b = 1. n = 10 ^, and m = 10 kms~^ to Fig. 6 
of lChieze p.tM lll99f^ . We find a very good agreement, 
though we use a very different method of integration, that 
cooling processes have been strongly refreshed to account 
for dissociative shocks, and that our resolution is higher 
at both shocks and relaxation layer. The method used by 
IChieze et all l)l998f) to recover the time is hence validated. 

In particular, we confirm the long time scales to steady 
C-type shocks found by these authors, of the order of the 
ion crossing time scale, and their weak dependence on the 
magnetic field for < 5 < 1. We add here (see Table ^ 
that these time scales are roughly inversely proportional 
to the density, as expected from the cooling time scales, 
and that there is no significant change with respect to the 
velocity, except when dissociative and n on-dissociative ve- 
lociti es are compared (see also Fig. 4 of iLe Bourlot et 
I2OO2I for example ). The length scales of the shocks 



which actually dissociate molecules are longer than the 
non-dissociating ones by a factor of 5-10, due to the fact 
that II2 is a very efficient coolant. 

3.3. Conditions for steady CJ-type shocks 

When the magnetic field is lower than a critical value, 
the steady shock is composed of a magnetic precursor, 
followed by a J-shock behind which the two fluids re- 
coupl e, and the ga s relaxes towards its final post-shock 
state. IChieze et ah (jl99&) have already shown one such 
model for non-dissociative velocities in a diffuse medium 
(for riH — 25 cm~^, u = 10 kms~^ and B = 5 ^G). This 
work presents models for dissociative and partly ionising 
velocities, which were out of reach of their anamorpho- 
sis method, because it could not resolve both the J-front 
and the H2 reformation zone. The final steady structure of 
such a shock can be seen in Fig. 2. It mixes all the char- 
acteristics of a C-type precursor with those of a partly 
ionising J-type shock (see Fig. 5). 

iLe Bourlot et oj] l)2002|) computed the critical values of 
the magnetic field and velocity for the transition from C to 
CJ-type behaviour. Their criterion to determine when a J 
front would remain in a stationary shock with a magnetic 
field is the presence of a sonic point in the neutral veloc- 
ity profile (in the shock frame). Such a sonic point is not 
encountered in any of our CJ-type models. Hence, the cri- 
terion previously used for CJ transition is too strong, and 
the range of parameters where stationary C-type shocks 
exist should be reduced compared to their results. 

In Paper II, we give a means of determining the fate 
(C or CJ) of a magnetic shock with the only help of a 
steady-state code for non-dissociative velocities. 

3.4. Chemically driven oscillations 

We have identified two (non mutually exclusive) situa- 
tions where oscillations occur : weakly dissociative shocks 
which undergo large rebounds driven by II2 reformation, 
and partly ionising shocks which show smaller oscillations 
probably driven by H ionisation. 

3.4.1. Oscillations due to H2 dissociation/reformation 

We find a narrow range of velocities, close to the dissoci- 
ation limit Ud, where shocks cycle between an expansion 
phase where the shock is dissociative, and a retreat phase 
where the shock is non-dissociative (an example is shown 
in Fig. 4). Even if u < Ud, it can happen that the initial 
entrance velocity in the shock u° — u + v is greater than 
Ud^. We define such shocks as weakly dissociative shocks. 
Fig. shows the first period of the shock with parameters 
b ~ 0, n ~ 10^ cm~^ and u ~ 25 kms~^. 

^ Because the compression factor in the atomic plateau af- 
ter a fully dissociative shock is around 6 (see Paper II), this 
happens for u > |iid 
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In the dissociative expansion phase, the front is fol- 
lowed by an isothermal atomic plateau at T ~10'* K. The 
length of the plateau is usually governed by O cooling fol- 
lowed by H2 reformation. When H2 reaches a high enough 
abundance, it becomes again the main coolant. At this 
point, the higher the H2 abundance, the faster it cools 
and is compressed, and the faster it is reformed. This run- 
away process yields a very sharp exit from the plateau (see 
Fig. 4a). 

When the compression factor at the end of the plateau 
increases, gets closer to u. Therefore, u° < Ud, and 
H2 survives the shock. Fast cooling ensues, and the shock 
collapses until it rebounds near the piston, with a total 
thickness on the order of the H2 cooling length. The cycle 
is repeated and the shock undergoes large oscillations (see 
Fig. 4d). Here, the H2 reformation instability proceeds ex- 
ac tly like the therma l instability of the 150 kms^^ model 
of iGaetz et all lll988tl or the type C (as defined by them) 
models of IWalder fc Folinil l(l99fil) : the gas condenses iso- 
barically in a shell as shown on Fig. 

To assess whether the bounces ("arches") are persis- 
tent or damped in the 6 = case, we computed at least ten 
periods with the help of a reduced network of 8 species (see 
Fig. 4d). The characteristic periods (in years) and typical 
amplitudes (in pc) of the arches (A) are given in tabled 
These time and length scales are found to depend strongly 
on the parameters n and u. They also depend on the con- 
trol parameters of the grid, as well as on properties of the 
cooling functions (like the coUisional excitation of CO by 
H), although the qualitative features remain unchanged. 

For weakly dissociative shocks with non-zero 6, the J- 
type feature initially behaves like the same shock without 
magnetic field, but the magnetic precursor rapidly decel- 
erates the neutrals velocity below the threshold for disso- 
ciation, and further rebounds are switched off. Hence we 
observe only one arch when a magnetic field is present. 
The period (in years) and amplitude (in pc) of this first 
arch (A') are given in tabled They are very similar to 
their values for 6 = 0. 

3.4.2. Oscillations due to H ionisation 

For higher velocities, temperatures are sufficient to start 
to collisionally ionise H in the plateau. As a result, the 
electron fraction increases progressively. If the ionisation 
fraction is high enough (i.e. if u*^ is large enough) Lyman 
cooling becomes dominant over Oxygen, and the end of 
the plateau is first triggered by this coohng. Once the 
temperature gets below 10* K, Lyman cooling drops and 
the gas enters a second plateau. The end of the second 
plateau is triggered by run-away H2 reformation as previ- 
ously described. An example of the thermal, chemical, and 
cooling structure of such a shock is shown in Fig. 5. We 
define shocks that present two such plateaux as "partly 
ionising shocks" . 

In these shocks, the length of the first plateau is gov- 
erned by the ionisation length of H. Since the temperature 
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Fig. 6. Temperature and pressure snapshots for the first 
period of the bouncing shock with parameters 6 = 0, 
n — 10^ cm^"^, and u — 25 kms~^. The plotting scale 
for the temperature is logarithmic while it is linear for the 
pressure. 

Tp inside the first plateau is quadratic in u (see equation 
29 of Paper II), and the ionisation length is exponential 
in 1/Tp, there is a very strong dependence of this length 
on the entrance velocity. This strong dependence seems to 
drive an oscillatory behaviour (denoted "o" in Table 1), 
with an amplitude on the order of the ionisation length. 
It appears to be damped over a few periods in most cases. 
Figure 5d illustrates this temporal behaviour. 

Partly ionising shocks are not exclusive of weakly dis- 
sociative shocks, and a few shocks present both character- 
istics (those marked with both an " A" and an " o" in table 

ID- 

The relaxation layer behind a partly ionising J-front 
with 6 7^ is very similar to pure J-type shocks (compare 
Fig. 2 and Fig. 5). The only difference is the onset of a typ- 
ical density where thermal pressure is relayed by magnetic 
pressure. At this point, and for sufficiently late times, the 
compression stops. This makes the length scales of the H2 
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reformation layer behind this point much greater, but the 
thermal and chemical evolutions are quite the same, as 
well as the oscillation properties. 

3.4.3. Previously known oscillations of shocks with 
chemistry 

iLim et all ()2002(l and ISmith fc RosctI |)2003|) also en- 
counter oscillations in their simulations of time-dependent 
shocks with chemistry. We compare our findings to their 
work; 

iLim et all tOO^ use explicit adaptive mesh refine- 
ment with chemistry splitted from hydrodynamics. They 
also follow H ionisation and H2 dissociation in a time- 
dependent way. However, they do not report any insta- 
bility in the n = 10'* and m = 30 kms~^ case, while 
we see strong bouncing oscillations. This difference could 
arise because they do not model H2 reformation on the 
grains, but only the reformation via the H~ process, which 
is much slower. They do identify various oscillating be- 
haviours in the course of their accelerated low-density 
shock, including oscillations of small amplitude and short 
time scales related to the ionisation of H when u ~ 50 
km s~*, which may be similar to ours. However, their re- 
sults are difficult to compare to our constant inflow speed, 
m uch higher density simu lat ions . 

RosenI l|200,^) present shock simulations with 
the same boundary conditions and densities as ours. Their 
inflow velocities are too high {u > 40 km s to ob- 
serve the dissociation/reformation bouncing oscillation, 
since their shocks are always dissociative. However, they 
observe a "quasi-periodic or chaotic collapse" driven by 
CO reformation for a very wide range of velocities (u from 
40 to 60 kms~*). The main differences with our work are 
their approximation of chemical equilibrium and optically 
thin emission for CO which makes CO cooling dominant 
in the atomic plateau, along with their coarser time-step 
control and lower resolution. 



4. Discussion 

We emphasise here a few points that remain to be inves- 
tigated. 

The opacity of the CO, H2O, and OH molecular cool- 
ings could be more detailed, including an LVG treatment 
in the rapidly compressing parts of the shocks. This would 
increase their importance compared to H2, though per- 
haps not enough to trigger the large instabilities found 
in the dissociative shock simulations of Smith fc RosenI 
1)2003(1 . Even at our extremely high resolution, we have 
also observed an effect of the grid parameters on the ampli- 
tude and period of oscillations. A linear analysis remains 
to be carried out to determine accurate periods without 
the caveats of the numerics. 

Next, the reformation of molecules in the relaxation 
layer showed the necessity of a treatment of the chemi- 
cal diffusion processes. We did not check the influence of 
the diffusion model on the reformation zone. Taking the 



diffusion of heat into account might change the thermal 
behaviour of our shocks. At higher velocities, the Lyman 
cooling and probably a few more atomic coolings will not 
be optically thin anymore and should also be treated more 
carefulh^^ 

As iLe Bourlot et all (|2002() pointed out with steady 
models, the time-dependent tracking of the populations 
of H2 has a large impact on the dissociation and hence 
the dynamical behaviour of shocks. This effect deserves 
more investigation in a fully time-dependent context. 

Oblique magnetic fields and grain dynamics have been 
completely neglected in this study. The lack of a veloc- 
ity for the grains may prove to be the main caveat of 



steady-state calculations suggest ( 


Ciolek & Robergfl2002l 


iFlower & Pineau des Foretdl2003 


). In a steady-sate con- 



portance of oblique magnetic fields, although it has not 
yet been investigated in any multifluid, time-dependent 
studies. However, our numerical technique should prove 
useful to model both these effects. 

5. Conclusions 

The adaptive moving grid technique proves to be a 
unique tool to model time-dependent magnetised molecu- 
lar shocks. It is currently the only algorithm able to model 
time-dependent dissociative shocks in presence of a mag- 
netic field. With th is tool, we validate previous results by 
IChieze et all l)l998l) and we investigate the formation of 
shocks in the conditions of protostellar j ets. 

We find in agreement with Pi neau d es Foret s et all 
(I1997I) : ISmith fc Mac Lowl IiQQtIi : IChieze et all lll998^ 
that C-shocks are steady after a very long time delay. 
We produce the first C-type and CJ-type shocks in the 
dissociative range of velocities. We point out that the oc- 
currence of a sonic point is not a valid criterion for the 
transition to CJ-type steady states. 

In our simulations, we also discover two oscillatory be- 
haviours, which are linked to H2 dissociation/reformation 
and to Hydrogen ionisation, respectively. The oscillations 
vanish when a strong magnetic field is present and af- 
ter a significant magnetic precursor has built up. The 
oscillation periods and amplitudes are found to depend 
strongly on the density and inflow speed. A detailed treat- 
ment of CO opacity is required to assess the reality of the 
CO-driven instabilitie s in dissociative shocks reported by 
ISmith fc RosenI ^l2003^ . 

In a companion paper (Paper II), we analyse our re- 
sults in comparison with truncated steady state models. 
We derive analytical relations as well as constructions for 
intermediate times of non-dissociative shocks with or with- 
out magnetic fields. 

Acknowledgements. We thank the referee (Pr. T.W. Hartquist) 
for his critical remarks on this paper, which lead us to clarify 
its scope and contribution to the field, and to greatly improve 
its presentation. 
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Fig. 1. Stationary structure at t = 2 10^ yr and shock trajectories of a J-type shock with n = 10^ cm~^, u = 20 kms~^, 
and 6 = 0. (a) Thermal structure : temperatures of the neutrals, ions and electrons. The crosses on the curve are the 
points where the flow time computed in the frame of the shock is equal to 10" yr, n being the number indicated 
near the cross. '+' signs stand for neutrals and 'x' for charges, (b) Cooling structure : magnitude of the main cooling 
processes, (c) Chemical structure : abundance relative to hydrogen nuclei of species of interest, (d) Shock trajectory : 
Position (thick curve) and velocity (thin curve) of the J-front against time, in the frame of the piston. Fig. 2. Same 
as Fig. 1 for a steady CJ-typc shock with n = 10^ cm""^, m = 40 kms~^, and h = 0.1, at t = 2 IQ^ yr. The dashed 
curve in (d) plots the position of the C-precursor against time. Fig. 3. Same as Fig. 2 for a steady C-type shock with 
n = 10^ cm-3, u = 20 kms-\ and 6 = 0.1, at t = 10^ yr. 
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Fig. 4 : weakly dissociative J-shock 



(4a) 



b=0, n=107cm ^ u=25 km/s 



Fig. 5 : partly ionising J-shock 




(5b) 
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Fig. 4. Same as Fig. 1 for a weakly dissociative J-shocli of parainct{;rs n = 10^ cm~^, u = 25 kms~^, 6 = at t = 50 yr. 
In panel (d), the thick solid line is the trajectory of the shock in a simulation with a reduced network of 8 species. The 
simulation with the whole network is plotted in a thick dashed line, up to the point where it stalled. Fig. 5. Same as 
Fig. 1 for a partly ionising J-shock of parameters n = 10^ cm~^, u = 40 kms~^, 6 = at t = 220 yr. 
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(J=10-3) 




(j=io-*) 




0=500 


(3 10-*) 


o=25 


(3 10-5) 


0=3 


(3 10-*^) 



6=1 n = 103 cm-3 n = 10* cm-3 n = 105 cm-3 
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Table 1. For each run, we provide the final steady-state of each shock (J, CJ, C) followed by the time scale in years 
at which steadiness is achieved. The total width of the shock at this point (in pc) is given in parentheses (in the case 
of CJ-type shocks, J is the length of the relaxation layer, and C is the length of the magnetic precursor). C> 700 
indicates that a J-front is still present when the code stalls a,t t = 700 yr. A and o are the average duration of one 
large arch or one small oscillation, when present. The associated length scales correspond to the typical amplitude of 
these oscillations. A' stands for a lone arch, o' stands for damped small oscillations. 



